# Replication File: Manuscript PSRM-OA-2015-0037
# Wojcik, Stefan: Why Legislative Networks? Analyzing Legislative Network Formation
# Oct 24, 2016
#install.packages("ergm")
#install.packages("latticeExtra")
#install.packages("stargazer")

# Please note that the number of observations will need to be added manually. 
#setwd("/Volumes/TINY CRYPT/papers/Working Projects/Dissertation/Main Work Folder/Analysis/PSRM/Replication Archive Final/Replication_Archive")
library(ergm)
library(stargazer)
library(coda)
library(latticeExtra)

# Start log:
sink("PSRM_replication_logERGM_final.txt")

load("PSRM_replication_image_networks_noIMP.RData")
#Main ERGM results: (TABLE 1)
#control.min = control.ergm(MCMC.burnin=10000, MCMC.interval=50, MCMC.samplesize=20000, MCMLE.maxit=30)
control.min = control.ergm(MCMC.burnin=1000000, MCMC.interval=200, MCMC.samplesize=300000, MCMLE.maxit=100) 
control.nobackfill = control.ergm(MCMC.burnin=900000, MCMC.interval=150, MCMC.samplesize=500000, MCMLE.maxit=100)

print("Main Results from Table 1:")

set.seed(12)
fit1 <- ergm(netcom.abs~edges+  absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MLE", control=control.min)

set.seed(12)
fit2 <- ergm(netsoc.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MLE", control=control.min)

set.seed(12)
fit3 <- ergm(netinfo.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MLE", control=control.min)

# Summarize the three models
stargazer(fit1, fit2, fit3, title = "TABLE 1", out="results_from_table_1.txt")
summary(fit1)
summary(fit2)
summary(fit3)

## INTERPRETATION OF EFFECTS:
geteffects <- function(model, control, treatment){
  coefs <- model$coef
  if(length(control)==length(coefs) & length(treatment)==length(coefs)) {
    contr.pred <- exp(coefs %*% control)/(1+exp(coefs %*% control))
    treat.pred <- exp(coefs %*% treatment)/(1+exp(coefs%*% treatment))
    effect.out <-  contr.pred-treat.pred
    return(list(control.prediction=contr.pred, treatment.prediction=treat.pred, effect=effect.out))
  }
  else {print("lengths do not match")}
}

#Effect of party assuming they work on the same floor and they share two mutual friends
control =   c(1, 11, 1, 1, 1, 0, 0, 0, 2) # not same party
treatment = c(1, 11, 1, 1, 1, 0, 0, 1, 2) # same party
print("PARTY EFFECTS: Control case is different party, treatment is same party")
print("Effects interpreted as: control condition probability is ___, treatment condition probability is ___. Control condition is ___ more/less likely than treatment condition.") 
print("Estimated probabilities of edges given same and different party: model 1 (table 1), discussed on p. 19")
geteffects(fit1, control, treatment) 
print("Estimated probabilities of edges given same and different party: model 2 (table 1), discussed on p. 19")
geteffects(fit2, control, treatment) 
#control prob is ___, which is ___ more/less likely than treatment
print("Estimated probabilities of edges given same and different party: model 3 (table 1, discussed on p. 19)")
geteffects(fit3, control, treatment) 
# control prob is ___, which is ___ more/less likely than treatment

#Age
control <-   c(1,  0, 1, 1, 1, 0, 0, 0, 2) #same age
treatment <- c(1, 11, 1, 1, 1, 0, 0, 0, 2) #difference in age as 11 years
print("Effects interpreted as: control condition probability is ___, treatment condition probability is ___. Control condition is ___ more/less likely than treatment condtion.") 
print("Control case is same age, treatment is difference in age: 11 years")
print("Estimated probabilities of edges given same and different age: model 1 (table 1), discussed on p. 19 and 20")
geteffects(fit1, control, treatment)
print("Estimated probabilities of edges given same and different age: model 2 (table 1), discussed on p. 19 and 20")
geteffects(fit2, control, treatment)
print("Estimated probabilities of edges given same and different age: model 3 (table 1), discussed on p. 19 and 20")
geteffects(fit3, control, treatment)

#Floor
print("Exponentiated coefficients of leg. floor: model 1 (table 1), discussed on p. 20")
exp(fit1$coef['nodematch.floor']) # times as likely to form ties
print("Exponentiated coefficients of leg. floor: model 2 (table 1), discussed on p. 20")
exp(fit2$coef['nodematch.floor']) # times as likely
print("Exponentiated coefficients of leg. floor: model 3 (table 1), discussed on p. 20")
exp(fit3$coef['nodematch.floor']) # times as likely
control =   c(1, 11, 1, 0, 1, 0, 0, 0, 2) # both non-matching floor
treatment = c(1, 11, 1, 1, 1, 0, 0, 0, 2)
print("Effects interpreted as: control condition probability is ___, treatment condition probability is ___. Control condition is ___ more/less likely than treatment condtion.") 
print("Control case is different floor, treatment is matching floor")
print("Estimated probabilities of edges given same and different floor: model 1 (table 1), discussed on p. 20")
geteffects(fit1, control, treatment) 
print("Estimated probabilities of edges given same and different floor: model 2 (table 1), discussed on p. 20")
geteffects(fit2, control, treatment) 
print("Estimated probabilities of edges given same and different floor: model 3 (table 1), discussed on p. 20")
geteffects(fit3, control, treatment) 

#Leadership - ::the last number was updated::
control <-   c(1, 11, 1, 1, 0, 1, 0, 0, 2) #both leadership
treatment <- c(1, 11, 1, 1, 0, 0, 0, 0, 2) 
print("Effects interpreted as: control condition probability is ___, treatment condition probability is ___. Control condition is ___ more/less likely than treatment condtion.") 
print("Control is both leadership, treatment is leadership to rank and file")
print("Estimated probabilities of edges given both leadership and mixed leadership/non-leadership: model 1 (table 1), discussed on p. 20")
geteffects(fit1, control, treatment)
print("Exponentiated coefficients of leadership: model 1 (table 1), discussed on p. 20")
exp(fit1$coef['nodematch.leadership.TRUE']) # 70% greater
print("Estimated probabilities of edges given both leadership and mixed leadership/non-leadership: model 2 (table 1), discussed on p. 20")
geteffects(fit2, control, treatment)
print("Exponentiated coefficients of leadership: model 2 (table 1), discussed on p. 20")
exp(fit2$coef['nodematch.leadership.TRUE']) # 70% greater
print("Estimated probabilities of edges given both leadership and mixed leadership/non-leadership: model 3 (table 1), discussed on p. 20")
geteffects(fit3, control, treatment) # has been updated
print("Exponentiated coefficients of leadership: model 3 (table 1), discussed on p. 20")
exp(fit3$coef['nodematch.leadership.TRUE']) # 70% greater

#State 
control=  c(1, 11, 1, 1, 1, 0, 0, 0, 2) # UF is state
treatment=c(1, 11, 1, 1, 1, 0, 1, 0, 2)
print("STATE EFFECTS: Control is different state, treatment is same state")
print("Effects interpreted as: control condition probability is ___, treatment condition probability is ___. Control condition is ___ more/less likely than treatment condtion.") 
exp(fit1$coef['nodematch.UF'])
print("Estimated probabilities of edges given same and different state: model 1 (table 1)")
geteffects(fit1, control, treatment) 
exp(fit2$coef['nodematch.UF'])
print("Estimated probabilities of edges given same and different state: model 2 (table 1)")
geteffects(fit2, control, treatment) 
exp(fit3$coef['nodematch.UF'])
print("Estimated probabilities of edges given same and different state: model 3 (table 1)")
geteffects(fit3, control, treatment) 

### FIGURES 7-15::: can be viewed in Rstudio. :::Too large to be uploaded to Dataverse:::

# PLOTTING PARAMETERS:
change_params = function(){
  oparm = par()$mar
  oparo = par()$oma
  par(mar=c(2, 2, 2, 2))
  par(oma = c(.2, .2, .2, .2))
  print("FYI: your margins were too small, attempting to change them. Figures may not look quite right.")
}

# MCMC diagnostics
#pdf("fit1.pdf")
mcmc.diagnostics(fit1)
#dev.off()

#pdf("fit2.pdf")
mcmc.diagnostics(fit2)
#dev.off()

#pdf("fit3.pdf")
mcmc.diagnostics(fit3)
#dev.off()

### FIGURES 16-27::: can be viewed in Rstudio. :::Too large to be uploaded to Dataverse:::

# GOODNESS OF FIT PLOTS
#pdf("gof_fit1_t1.pdf")
tryfit1 = try(plot(gof(fit1)), silent = T)
if(class(tryfit1) == "try-error"){
  change_params()
  tryfig1 = try(plot(gof(fit1)), silent = T)
}
#dev.off()

#pdf("gof_fit2_t1.pdf")
tryfit2 = try(plot(gof(fit2)), silent = T)
if(class(tryfit2) == "try-error"){
  change_params()
  tryfig2 = try(plot(gof(fit2)), silent = T)
}

#dev.off()

#pdf("gof_fit3_t1.pdf")
tryfit3 = try(plot(gof(fit3)), silent = T)
if(class(tryfit3) == "try-error"){
  change_params()
  tryfig3 = try(plot(gof(fit3)), silent = T)
}
#dev.off()

# MPLE Results for Comparison (TABLE 5)
set.seed(12)
print("Main Results from Table 5:")
fit1 <- ergm(netcom.abs~edges+  absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MPLE")
fit2 <- ergm(netsoc.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MPLE")
fit3 <- ergm(netinfo.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MPLE")
stargazer(fit1, fit2, fit3, title = "TABLE 5", out="results_from_table_5.txt")
summary(fit1)
summary(fit2)
summary(fit3)

# GWESP TERMS: (Reviewer Materials Only)
print("GWESP TERMS: (Reviewer Materials Only). These do not appear in the paper.")
set.seed(12)
fit1 <- ergm(netcom.abs~edges+  absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + gwesp(0.2, fixed=T), estimate="MLE")
set.seed(12)
fit2 <- ergm(netsoc.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + gwesp(0.2, fixed=T), estimate="MLE")
set.seed(12)
fit3 <- ergm(netinfo.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + gwesp(0.2, fixed=T) , estimate="MLE")
stargazer(fit1, fit2, fit3, title = "REVIEWER MATERIALS ONLY", out="GWESP_rev_only.txt")

# Alternative Specification (TABLE 2)
print("Alternative Specifications: (Table 2)")
netcom.abs <- delete.vertices(netcom.abs, which(is.na(netcom.abs %v% "ideology")))
set.seed(12)
fit1b <- ergm(netcom.abs~edges+ absdiff("n.mandates") + absdiff("ideology") + nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F)+ transitiveties(), estimate="MLE", control=control.min)
netsoc.abs <- delete.vertices(netsoc.abs, which(is.na(netsoc.abs %v% "ideology")))
fit2b <- ergm(netsoc.abs~edges+ absdiff("n.mandates") + absdiff("ideology") + nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F)+ transitiveties(), estimate="MLE", control=control.min)
netinfo.abs <- delete.vertices(netinfo.abs, which(is.na(netinfo.abs %v% "ideology")))
fit3b <- ergm(netinfo.abs~edges+ absdiff("n.mandates") + absdiff("ideology") + nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F)+ transitiveties(), estimate="MLE", control=control.min)
stargazer(fit1b, fit2b, fit3b, title = "TABLE 2", out="results_from_table_2.txt")
summary(fit1b)
summary(fit2b)
summary(fit3b)

rm(list=ls())


# NO BACK-FILLING MODELS (TABLE 6)
print("Models that do not back-fill: (TABLE 6)")
load("PSRM_replication_image_networks_NOBACKFILL.RData")
control.nobackfill = control.ergm(MCMC.burnin=900000, MCMC.interval=150, MCMC.samplesize=500000, MCMLE.maxit=100)
set.seed(12)
fit1 <- ergm(netcom.abs~edges+  absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MLE", control=control.nobackfill)
fit2 <- ergm(netsoc.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MLE", control=control.nobackfill)
fit3 <- ergm(netinfo.abs~edges+ absdiff("age") +nodematch("education", diff=F)+nodematch("floor", diff=F) + nodematch("leadership", diff=T) + nodematch("UF", diff=F) +nodematch("party", diff=F) + transitiveties() , estimate="MLE", control=control.nobackfill)
stargazer(fit1, fit2, fit3, title = "TABLE 6", out="results_from_table_6.txt")
summary(fit1)
summary(fit2)
summary(fit3)

# Selection model: table 7 
# Making the selection bias table
#selection = read.csv("PSRM_selection_bias.csv", sep = "\t")
#mod <- glm(responded~., family = binomial(link="logit"), data=selection)
#stargazer(mod, title = "Selection Bias Model", out="results_from_table_7.txt")

sink()
